Chapter 7 Beta diversity
7.1 Captive-bred vs wild
wild_pre_samples <- sample_metadata %>%
filter(diet!="Post_grass") %>%
dplyr::select(sample) %>% pull()
genome_counts<- genome_counts %>%
column_to_rownames("genome") %>%
select(all_of(wild_pre_samples))%>%
rownames_to_column("genome")
sample_metadata <- sample_metadata %>%
filter(diet!="Post_grass")beta_q0n <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
genome_counts <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)
beta_q1p <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
remove_rownames() %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
beta_q1f <- genome_counts_filt %>%
# remove_rownames() %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)7.1.1 Permanova analysis
Richness
sample_metadata_row<- column_to_rownames(sample_metadata, "sample")
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]
betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.004224 0.0042240 7.8055 999 0.003 **
Residuals 22 0.011905 0.0005412
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Wild
Pre_grass 0.008
Wild 0.010583
adonis2(beta_q0n$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.01526961 | 0.2489907 | 7.29391 | 0.001 |
| Residual | 22 | 0.04605641 | 0.7510093 | NA | NA |
| Total | 23 | 0.06132602 | 1.0000000 | NA | NA |
Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.006925 0.0069254 0.5404 999 0.491
Residuals 22 0.281957 0.0128162
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Wild
Pre_grass 0.476
Wild 0.47005
adonis2(beta_q1n$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 1.538238 | 0.2539205 | 7.487472 | 0.001 |
| Residual | 22 | 4.519715 | 0.7460795 | NA | NA |
| Total | 23 | 6.057953 | 1.0000000 | NA | NA |
Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.007671 0.0076708 1.8217 999 0.179
Residuals 22 0.092637 0.0042108
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Wild
Pre_grass 0.187
Wild 0.19084
adonis2(beta_q1p$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.3720264 | 0.3082818 | 9.804858 | 0.001 |
| Residual | 22 | 0.8347476 | 0.6917182 | NA | NA |
| Total | 23 | 1.2067741 | 1.0000000 | NA | NA |
Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.003070 0.0030698 1.7361 999 0.171
Residuals 22 0.038901 0.0017682
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Wild
Pre_grass 0.176
Wild 0.20119
adonis2(beta_q1f$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.01079711 | 0.187686 | 5.083124 | 0.103 |
| Residual | 22 | 0.04673041 | 0.812314 | NA | NA |
| Total | 23 | 0.05752752 | 1.000000 | NA | NA |
7.1.2 Beta diversity plots
7.1.2.1 Richness
beta_q0n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet", "Wild" = "Wild"))+
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(face = "bold", size = 12),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.position = "right", legend.box = "vertical"
) +labs(color='Origin')7.1.2.3 Phylogenetic
beta_q1p$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Origin')7.1.2.4 Functional
beta_q1f$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Origin')7.2 Before and after grass
genome_counts <- genome_counts_filt
post_pre_samples <- sample_metadata %>%
filter(diet!="Wild") %>%
dplyr::select(sample) %>% pull()
genome_counts<- genome_counts %>%
column_to_rownames("genome") %>%
select(all_of(post_pre_samples))%>%
rownames_to_column("genome")
sample_metadata <- sample_metadata %>%
filter(diet!="Wild")beta_q0n <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
genome_counts <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)
beta_q1p <- genome_counts %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
remove_rownames() %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)7.2.1 Permanova analysis
Richness
sample_metadata_row<- column_to_rownames(sample_metadata, "sample")
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]
betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0000047 0.00000467 0.0151 999 0.888
Residuals 22 0.0067978 0.00030899
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Post_grass
Pre_grass 0.885
Post_grass 0.90328
adonis2(beta_q0n$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.0006002558 | 0.030093 | 0.6825872 | 0.785 |
| Residual | 22 | 0.0193464336 | 0.969907 | NA | NA |
| Total | 23 | 0.0199466894 | 1.000000 | NA | NA |
Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.001983 0.0019829 0.1447 999 0.686
Residuals 22 0.301550 0.0137068
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Post_grass
Pre_grass 0.686
Post_grass 0.70733
adonis2(beta_q1n$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.147998 | 0.03074149 | 0.6977632 | 0.946 |
| Residual | 22 | 4.666278 | 0.96925851 | NA | NA |
| Total | 23 | 4.814276 | 1.00000000 | NA | NA |
Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.001772 0.0017715 0.4816 999 0.468
Residuals 22 0.080927 0.0036785
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Post_grass
Pre_grass 0.465
Post_grass 0.49497
adonis2(beta_q1p$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.02764442 | 0.04145548 | 0.951464 | 0.59 |
| Residual | 22 | 0.63920164 | 0.95854452 | NA | NA |
| Total | 23 | 0.66684606 | 1.00000000 | NA | NA |
Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.008352 0.0083519 3.722 999 0.037 *
Residuals 22 0.049367 0.0022439
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Pre_grass Post_grass
Pre_grass 0.03
Post_grass 0.066699
adonis2(beta_q1f$S ~ diet,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| diet | 1 | 0.002870065 | 0.03722968 | 0.8507252 | 0.452 |
| Residual | 22 | 0.074220722 | 0.96277032 | NA | NA |
| Total | 23 | 0.077090787 | 1.00000000 | NA | NA |
7.2.2 Beta diversity plots
7.2.2.1 Richness
beta_q0n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Origin')7.2.2.3 Phylogenetic
beta_q1p$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Origin')7.2.2.4 Functional
beta_q1f$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(diet) %>%
mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
geom_point(size = 4) +
scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
) +
labs(color='Origin')